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APPLICATION OF A ONE-STRIP INTEGRAL METHOD TO THE UNSTEADY 
SUPERSONIC AERODYNAMICS OF AN INCLINED FLAT SURFACE* 

By Robert M. Bennett 
Langley Research Center 

SUMMARY 

The method of integral relations is applied in a one-strip approximation to the per- 
turbation equations governing small motions of an inclined, sharp-edged, flat surface 
about the mean supersonic steady flow. Algebraic expressions for low- reduced-frequency 
aerodynamics are obtained and a set of ordinary differential equations are obtained for 
general oscillatory motion. Results are presented for low-reduced-frequency aerody- 
namics and for the variation of the unsteady forces with frequency. 

The method gives accurate results for the aerodynamic forces at low reduced fre- 
quency which are in good agreement with available experimental data. However, for cases 
in which the aerodynamic forces vary rapidly with frequency, the results are qualitatively 
correct, but of limited accuracy. Calculations indicate that for a range of inclination 
angles near shock detachment such that the flow in the shock layer is low supersonic, the 
aerodynamic forces vary rapidly both with inclination angle and with reduced frequency. 

INTRODUCTION 

Recent developments in flight at high Mach numbers and in computer technology 
have led to several highly developed methods for the analysis and computation of complex 
steady flow fields. For example, the finite difference methods (refs. 1 to 3) and the 
method of integral relations (refs. 4 to 7) have been applied to the analysis of complicated 
flows over blunt and pointed bodies. However, the surveys of references 8 to 13 indicate 
that the corresponding unsteady aerodynamics for analyzing flutter, loads, and stability 
have not been as well developed. 

Linearized flow theory was used in much of the early analysis of unsteady super- 
sonic flows (refs. 14 and 15). Recently refined and implemented for an extended range 
of planforms (ref. 16), linearized theory is still inherently limited to low supersonic 
Mach numbers and small deflections. The development of piston theory (refs. 17 to 19) 
permitted extensions to higher Mach numbers, and this theory has been widely applied 
(refs. 20 and 21, for example). However, piston theory is based on hypersonic small- 
disturbance theory (ref. 22) arid is thus limited to small slopes. It is restricted to values 

*The material presented herein is based on a thesis entitled "An Analytical Investi- 
gation of an Oscillating Wedge in a Supersonic Perfect Gas Flow" submitted in partial ful- 
fillment of the requirements for the degree of Doctor of Philosophy in Aerospace Engi- 
neering, North Carolina State University at Raleigh, December 1971. 



of Mach number times deflection angle of less than one, and in the simplified version 
usually applied, it does not account for the details of the local flow. Other methods used 
include the simple Newtonian and shock- expansion theories (refs. 23 to 25), but they are 
also highly simplified and restricted in their range of applicability. Considerable 
research has been devoted to the refined analytical treatment of the unsteady aerody- 
namics of flat, sharp- edged surfaces, such as wedges or cones, with simple steady flow 
fields, in an effort to gain some insight into the limits of the simpler methods and into the 
effects of strong shock waves and wave interactions on the unsteady flows (refs. 26 to 41). 
Analysis of the more complex configurations encountered in practice by the methods used 
for these simpler cases would, however, appear to be prohibitively difficult. 

The purpose of this report is to examine the application to unsteady flows of one 
method that has been used for the computation of steady flow fields — the method of inte- 
gral relations (refs. 4 to 7) - and to examine some pertinent trends for the unsteady aero- 
dynamic forces. 

The method of integral relations is essentially a method for reducing a system of 
partial differential equations of two (or more) independent variables to an approximating 
system of ordinary differential equations in one independent variable. Herein, application 
of the method is made in a one- strip approximation to the unsteady aerodynamics of an 
inclined, sharp-edged surface such as the compression side of a wedge or flat plate. Both 
the previously mentioned analytical results and some low-frequency experimental data 
(ref. 42) for the oscillating wedge are available for comparison. The linearized perturba- 
tion equations describing infinitesimal oscillatory motions about the steady flow condition 
are treated. Linearization and treatment of a specific motion reduces the number of inde- 
pendent variables from three (two space variables plus time) to only the two space varia- 
bles with frequency as a parameter. The method of integral relations is then applied to 
further reduce the governing equations to ordinary differential equations in one independent 
variable. The analysis is outlined and sample results are presented and discussed, with 
emphasis on the effects of the approach to detachment, surface inclination angle, and ratio 
of specific heats y. The oblique-shock perturbation boundary conditions, the calculation 
of coefficients for a symmetrical wedge from the coefficients for a single surface, and a 
computer program listing are given as appendixes. 

SYMBOLS 


a local speed of sound 

B 0 = tV - 1 
C chord force 


2 



pressure coefficient, 


C p 

c 

D 

e 

hx 

h y 

i = fi 


p - p. 


de- 


reference length 

quantity defined by equation (10) 

base of natural system of logarithms, 2.71828 

displacement parallel to mean surface position (fig. 1) 

plunging displacement normal to mean surface position (fig. 1) 


reduced frequency, ■ 

2V to 

motion-related aerodynamic forces and moments measured about the leading 
edge for a single surface (see eq. (11a) and eq. (lib), respectively) 

L^,M^ motion- related aerodynamic forces and moments measured about the leading 
edge for a symmetrical wedge (see eqs. (Bl)) 

Lj,Mj coefficients in normal-force and moment expressions for a single surface, 

j = 1,2, 3, 4, 7, 8 


k 

L',M’ 


Lj )W ,Mj ;W coefficients in normal-force and moment expressions for a symmetrical 
wedge, j = 1,2, 3, 4, 7, 8 


M 


Mach number, — — 


P V 


8v n 


P 




9Pg 

»/3 x 


p pressure 

q dynamic pressure 


3 




R 

0 " 9 0i 


T 

T 0 


shock- layer thickness in y-direction 

shock- layer thickness in y-direction for steady flow 


t time 


U v = 


9Ug 

9v n 



8Ug 

9^1 


a velocity component in x-direction (fig. 1) 


9V, 

V = - 

V V 9 V . 


n 


free-stream velocity 


V 


/3 = 


!l6 

0j8l 


velocity component in y-direction (fig. 1) 


perturbation of velocity of shock wave measured normal to steady shock 
wave (fig. 1) 


x,y 


x o»y 0 


^0 


surface coordinates (fig. 1) 

pitch-axis coordinates (as fractions of c), with x Q measured from leading 
edge and y Q measured from chord 

shock wave angle for steady flow (fig. 1) 
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H 

y 

«1 

6w 

e 

0 

© d 

9 

A 0 

P 

T 

U> 

Subscripts: 

s 

u,Z 

w 

6 

0 


perturbation of shock wave angle (fig. 1) 
ratio of specific heats 

perturbation of shock layer thickness, measured normal to steady shock 
wave (fig. 1) 

wedge angle ^6 W = - tan' *r for symmetrical wedge airfoil j 
perturbation parameter 
surface inclination angle 

surface inclination angle for shock detachment 
pitch angle 

included angle between shock and surface for steady flow (fig. 1) 
density 

ratio of airfoil thickness to chord c 
circular frequency of oscillation 

pertaining to perturbation motions of h x , hy, or 9, respectively 
denotes complex amplitude of perturbation quantity evaluated at surface 
pertaining to upper and lower surface, respectively, for a symmetrical wedge 
pertaining to a symmetrical wedge 

denotes complex amplitude of perturbation quantity evaluated just aft of 
shock wave 

steady flow value in shock layer 
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1 


denotes complex amplitude of perturbation quantity 


°o free-stream value 

Superscript: 

' pertaining to pitch (or pitching moment) about leading edge 

A bar over a symbol indicates a dimensional quantity. A tilde (~) denotes, forces 
normal to the surface. 


ANALYSIS 

A two-dimensional flat surface with a sharp leading edge and exposed to a super- 
sonic free-stream flow of a perfect gas at a mean inclination angle of 0 is considered 
while undergoing a specified time -dependent motion (fig. 1). Three types of simple- 
harmonic oscillations are considered separately - pitching about the leading edge 9 , 
translation normal to the mean position of the surface hy, and translation parallel to the 
mean position of the surface h x . The motion is considered to be an infinitesimal per- 
turbation about the mean or steady-flow condition and leads to a linearized analysis 
describing the region between the moving surface and shock-wave boundaries beginning 
at the leading edge and extending a preassigned distance downstream. Linearization of 
the governing equations permits the application of the boundary conditions at the mean 
position of the moving boundaries and also permits treatment of any rigid-body, time- 
dependent motion of the flat surface by superposition of the three motions considered. 
There are essentially three criteria that govern the amplitude limits for small oscilla- 
tions: (1) fraction of = ^0 " (2) fraction of 0 (or for very small 0, fraction of 

p^), or (3) fraction of 0^-0 where 0^ is the detachment angle. It might be noted 
that only compression surfaces are treated; however, in a linearized analysis of an expan- 
sion surface, the linear theory of reference 14 can be applied by considering the mean 
surface flow to be the effective free stream. Also for application to wedges, the base 
pressure is assumed constant and thus does not enter into the motion-related force and 
moment coefficients. 


Development of Perturbation Equations 

The governing flow equations written in divergence form for the x,y coordinate sys- 
tem (for example, reference 43, chapter 7) are as follows: 
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Continuity: 


dp d(pu) 8(pv) 

— + ■ ■ ■"__ - + = 0 
8t 9x dy 


x-momentum: 

8 (pu) 8 (p + pu 2 ) 8 (puy) 

9t 9x 8y 

y -momentum: 

9(pv) 9(pu v) a(p + pv 2 ) 


Energy: 


— f- — 4* 

at 9x 


9y 


= 0 


_a_ 

at 


r + ^7^ p ^ u2 + v2 ) + ^f p + ^^-^ 2+v2 ) 


(la) 


(lb) 


(lc) 


+ #=<^ 

9y 


p + (y- 1 ) —i—2 . — 


2y 


p(u 2 + v 2 )\=0 


(Id) 


where the barred quantities are dimensional. 

Small variations about the steady flow are considered by substituting in equations (1) 

^7\ ( v" t 


P = P Q + ep 1 (x,y)e 1 

p = p Q + 

u = u 0 + eir 1 (jr,y)e lt,,t 


v = ev 1 (x,y)e la;t 


where e is a small quantity; aJ is the frequency of oscillation; p^, p^, uj, and v^ 
are the complex amplitudes of the perturbations in the flow variables; and p Q , p Q , and 
Uq pertain to the steady or mean wedge flow. Neglecting terms of order of g2 an( j 
higher gives, for the y-momentum equation, for example, 


9vi 9pi 

p o u o-ir + if +1Ct,p o v i = 0 
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The equations are simpler in form if normalized by the shock-layer variables Pq and 
u 0 and if the shock-layer wavelength u 0 /w is used as the characteristic length. By 
setting p x = Pi/(p 0 u 0 2 ), P x = JS\ /Pq> = "l/^O’ V 1 = x = xw/uq, and 

y = yZ5/viQ, the following equations result: 

y-momentum: 




0 


x- momentum: 


8uj 

IF 


+ 


f£l 

9x 


+ 1U 


1 


= 0 


Energy: 


9 Pi 
dx + 


dvi 


M, 


2 ay 


■+ l 


O' Pi ~ U 1 - Q 


Br 


Continuity: 


(2a) 


(2b) 


(2c) 


9 Pl 
ax + 




= 0 


(2d) 


where Bq 2 = Mg 2 - 1. 

Note that the equations have been manipulated to eliminate p^ from equations (2a) 
to (2c). Equation (2d) is needed only to calculate p^ and will be disregarded here. Thus 
equations (2a) to (2c) are a system of three simultaneous, linear partial differential equa- 
tions for Pp u^, and vj as functions of x and y and with complex, constant 
coefficients . 


Boundary Conditions 

Small variations about the steady flow have been considered by substituting in equa- 
tions (1) for p, for example, 

P = P 0 + ep^yje 1 ^ 

where pj is a complex amplitude and e is the real perturbation amplitude parameter 
B, h x , or hy. The factor ee^*, common to all terms and thus deleted in equations (2), 
also appears in all perturbation boundary conditions and is omitted so that the boundary 
conditions are specified in terms of the complex amplitudes. 
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II III 


IB 


I 


II 


Boundary conditions at the leading edge.- For the assumed attached leading- edge 
shock, the shock-displacement boundary condition is that the shock wave displaces with 
the leading edge. In terms of the complex perturbation amplitude, this boundary condition 
for pitching motion is 

6 X (0) = 0 

(3a) 

for plunging motion normal to surface, 


6i(0) = -cos 

(3b) 

and for translation parallel to surface, 


6-^(0) = sin Aq 

(3c) 

Boundary conditions at the surface.- The tangent flow condition at the surface deter- 
mines the complex amplitude at the surface. For infinitesimal pitching motion, the bound- 
ary condition is 

v s = -(1 + ix) 

(4a) 

for plunging motion normal to the surface, 


v s = 

(4b) 

and for translation parallel to the surface, 


v s = 0 

(4c) 

Boundary conditions at the shockwave.- The perturbations of 
aft of the moving oblique shock wave are discussed in appendix A. 
tudes of the flow variable perturbations are given by equations (A2) 

the flow variables just 
The complex ampli- 
in the following form: 

P 6 = *pPl + iP v 6 i 

(5a) 

u 6 = + iU v 5! 

(5b) 

v 6 = V ^1 + iV v 6 l 

(5c) 

where 


d6i 

^i = d^ cosX o 

(5d) 

and where = 8ug/9/3i, and so forth, are derivatives of the shock rela- 

tions evaluated at the steady or mean conditions. 
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Approximate Solution by the Method of Integral Relations 


Development o f one -strip equations.- For the one-strip integral-relations approach 
used herein, a linear y-variation of Pj, u x , and v x is assumed. For example, 


P x (x,y) = P s (x) + [p 6 (x) - P s (x)J^ 


( 6 ) 


where subscripts s and 6 denote complex perturbation amplitudes evaluated at the 

6*(x) 

surface and just aft of the shockwave, respectively, and T(x) = T n + — is the total 

u cos Xq 

shock-layer thickness. Equation (6) and similar relations for u^ and v x are used in 
equations (2a) to (2c) and integrated in the y-direction. Products of perturbation quanti- 
ties are neglected. Sample terms for p^ are 


f T 

j 0 Pi (x >y) d y = [p s < x > 


+ p 6 < x) ] 


T 0 W 

2 


£ 


T 3p x (x,y) 

— ^y - Fr 

o 9 y 0 


dy = p*(x) - p„(x) 


and, with the use of Leibniz’ rule, 


C ^T^ dy - slj Pi(x,y)dy - fflp 6 (x) 

ow ir / \ / M dT o( x ) 

— + |[p s ^ ' P 6 (x) ]-d^- 


dpg(x) 

T 0 (x) 

dx 

2 


The results of this procedure are as follows: 

arf. a + P s + u 6 + u s) - ^>6 - e s + “6 - + ‘(»6 + «s) - 0 

Sf'e + v s) - - v s)^T + - p s) + *( v 6 + T s) - 0 

d 1 dTn 2 

^( p 6 + p s ) - ^(p 6 - p s )-ar + - v s) 
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+ B^[ M ° 2 ( P6 + PfS ) ’ + U s)] " 0 



Now, 


Tq(x) = x tan Xq 


and 


dT Q (x) 

dx 


= tan Xq 


Furthermore, the terms subscripted by 6 are given by the boundary conditions just aft 
of the shock wave (eqs. (5a) to (5c)). For example, 


p 6 (x) = P| 3 /3 1 (x) + iP v 6 1 (x) 


and thus 


dpg 

dx 



d/3 x (x) 

dx 


+ iP V 



d/3j(x) 


dx 


+ l 


Py 

COS Xq 


J3 x (x) 


as = _JL cos Xq (eq. (5d)). Developing and using similar expressions for du g /dx 
and dvg/dx gives the following equations: 


d/3i 

x T^- + 
dx 


2P 


P 


cot Xq - 1 + ix 1 + 


Vp cos X 0/ 




+ i 


2P V , „ . . V v 

— cot X 0 - (1 - ix) — 




s 


Vg dx V„ 


x ^r + (1 + “K + x tr + p s + X ( P |3 + U /3)^T 


+ (1 - ix)U fi - ix 


Py + Uy 


P COS Xr 


^1 - i[Py + Uy(l - ix)]6 x = 0 


(7a) 


(7b) 
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L 


x 5-^-+ (l + ix ^Qp - i-^-u s + xP fi ^ 
dx l B 0 2 s B 0 2 S ^ dx 


2V 


0 


B 0 2 tan Xq 


P |1 . ix ?\) + 

0\ R 2/ cos X, 
B 0 




o B 0 


h 


+ 1 


2V t 


Mr 


B n z tan X/- 


P vl X " te ~ t) “ “T 

1 B 0 / B o 2 


6 2v s 

1 BQ2tan Xq 


(7c) 


d6j Pi 


dx COS Xq 


(7d) 


The partial differential equations (eqs. (2)) and their boundary conditions have been trans- 
formed to four inhomogeneous, complex, variable -coefficient, ordinary differential equa- 
tions in p g , u s , 6 j_, and j 3 j. The known surface velocity v s (x) appears as a forcing 
function, and the required initial conditions are also contained in equations (7) as will be 
subsequently shown. 

Initial values at the leading edge.- Taking the limit as x - 0 in equations (7) gives 


ei<°> ■ ^ 


- 1 r- M°> 

0 V 0 


(8a) 


Pg(0) = P^O) + iPyS^O) 

Ug(0) = 11^(0) + iUySjCO) 


(8b) 

(8c) 


These conditions are the exact quasi-static linearized conditions from oblique-shock 
theory. 

Initial derivatives at the leading edge.- Differentiating equations (7) with respect to 
x, taking the limit as x - 0, and applying equations (8) gives an algebraic system which 

d/3 t dp du 

can be solved for — - , -r-s- , and 

x=0 x=0 x=0 
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The results are as follows: 


d£i 

dx 


Bg^tan Xq 


v R (0) 




+ i( M 0 2 Py - U v + V v ^ tan X 0 )«j(0) + ( M^P ' p - V 0 + 


+ -s-Vo tan X n + 


W 




(9a) 


dp s 

dx 


= P 


x=0 


d/3i 
0 dx 


+ tan X, 


x=0 


dVe 
0 dx 


+ 1 


tan Xq 


|x=0 


v s (0) + iVyfijW 


and 


- ( v f - ^)% (0 > 


(9b) 


dUg 

dx 


lx=0 


.^Es 

dx 


lx=0 


iu s (0) 


(9c) 


In equation 9(a), 


D = V|3 + p /3 B 0 2tan X Q 


( 10 ) 


Definition of force and moment coefficients .- The motion- related aerodynamic force 
and moment coefficients are defined for a single (upper; see fig. 1) surface by 


l' - - C [p(*) - P 0 ]<R- = 45-„iTk2 2 (l, + ltj)!* 


+ (Lg + i£j)* + 2^.7 + iLgJ^ 


,iwt 


(11a) 
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(lib) 



In equations (11), moments and displacements are measured at the leading edge. These 
definitions are similar to those often used in two-dimensional supersonic flutter analysis 
(e.g., ref. 14). Here, however, the force is normal to the surface (indicated by ~ over 
the L), rather than in the lift direction, and subscripts 7 and 8 pertain to fore and aft 
translation parallel to the mean surface position. (The symbol N and the subscripts 5 
and 6 are often used for unsteady flap hinge moments and for pertaining to a flap, respec- 
tively; see ref. 14.) It might also be noted that no forces along the surface are generated, 
since the effects of viscous skin friction have been neglected and the base pressure has 
been assumed to remain constant. In appendix B the aerodynamic coefficients for a wedge 
are given in terms of these definitions (eqs. (11)). 

Low-k solution b a sed on initial con ditions.- The independent variable x ^that is, 
x o^u serves as a nondimensional length scale for the oscillating wedge. For a specific 
value of x (such as c or unity), k = ^Tq^V^x and thus small values of x imply 
small values of k. The surface pressure for small values of x (and hence k) can be 
approximated as 


P s (x) = P s (0) + ^-1 


x=0 


( 12 ) 


Using the boundary conditions given by equations (3) in equations (8) and (9) and integrating 
for the force and moment coefficients as defined by equations (11) gives for small values 
of k and for a single exposed surface 


• m; 


P0 d Ps,hy 

i > 

2^oo dx 

x=0 


p o u o 


kL 9 = kM« s - _ p h (0) 

2 2 2p JT F s,h y 


k2L' = k2Mg = - — ° U ° p q (0) 

3 3 2^2 B,9 


(13) 


(Equations continued oa next page) 
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II I I II I II I II !■ 1111 I 


i lain mi ■■ 




p 0 *0 dp s , 0 


2 Pco U co ^ 


lx=0 


f - 3 M > p 0 dPs > h x 
L7- 4 M 7- _ 2P dF^ 

OO 


x=0 


~ Pn u n 

kL 8 = ™8“-^E's,h x ( ( » 


Evaluating equations (13) using equations (3), (4), (8), and (9) leads to the following alge- 
braic expressions for the coefficients involving only the properties of the steady flow p^, 
Uq, and so forth, and the shock-wave derivatives P^, V v , and so forth: 


where 


kL„ = — ^ 0 _ 0 — /p _ A cog x \ 

2 2puV fl U 1 0 

^ OO OO p 


k2t , . >W, 2 . p 


Pn u n 

kLo = — Ai sin X n 

8 VJUVf 1 0 


1, = ° | (V v + M„2p. sin .UAj - PpAj sin X 0 1 

P oo p* *- ^ 

Lj = -cot A 0 L 7 - ^ 


p 0 A 3 


kL’ = 


PQ u O 


4 2p u 

• OO 


OO oo 


tan Xq + 


P^(l - B 0 2 tan 2 Xo) 


D 


+ Ac 


Al = - P v V^g 

A 2 = U^Vv - U v V^ 


(14) 
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f -I Ai 
A3 = tan A 0 <1 - — i— — — - 
0 u> \ VoD sin A 

and D is given by equation (10). 

Numerical solution of one- strip equations.- To determine the variation of the force 
and moment coefficients with k, equations (7) are numerically integrated by using a com- 
puter subroutine for first-order, real, ordinary differential equations. These complex 
equations (eqs. (7)) are converted to a system of eight coupled, first-order equations 
relating the real (Re) and imaginary (Im) parts of 6j, j 8 j, u g , and p g . Initial values 
and derivatives at a small value of x are calculated from equation (12), for example, 
which requires using equations (8) and (9) for starting the numerical integration of the 

UqX 

differential equations. As k = -~- = ■■ _ , running integrals over x of the real and 

2V ao 2Voo 

imaginary parts of the pressure give the variation of the aerodynamic coefficients with 
k. By using the definitions given as equations (11), the coefficients are related to the 
integral of the pressure for a single exposed surface as follows: 
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where £ is the variable of integration in the x-direction. Similarly, for pitch, 

P f ,T - 2 

2p 


k2£ 3 = 'i?li.C Re M d « 


oo 00 



k 2 M ' 3 = 



4 d^ 


kM^ = 


2 Pq u q 

PooV 3 





A computing program is given in appendix C which calculates twice the above integrals. 

The trapezoidal rule is used for numerical integration of the surface pressures in the 
integrals to determine the coefficients. 

RESULTS AND DISCUSSION 
Comments on the Method 

One result of applying the method of integral relations is a change in the region of 
influence of the governing differential equations. Disturbances propagate along the normal 
coordinate rather than along characteristics and the description of wave-type phenomena 
is thus altered. The application to steady supersonic flow over pointed bodies has been 
discussed in reference 7. It was found that although the results of the method differed in 
detail, the exact wave behavior was approximated as a result of compensating effects. 
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The integral method has also been applied to the perturbation of plane entropy layers 
which involves extensive wave effects with favorable results (ref. 6). As discussed in 
references 29 and 37, the acoustic waves generated by the surface and shock wave are 
important factors in determining the unsteady forces on the wedge. The description of 
these wave phenomena by the one-strip integral method (eqs. (7)) is investigated by apply- 
ing the method to a steady angular perturbation of a wedge surface. The results are com- 
pared in figure 2 with the exact solution from hypersonic small-disturbance theory of 
reference 29. The results of the one-strip method do not follow the exact wave pattern, 
but approximate it in a smoothed oscillatory fashion. This is thus an important limita- 
tion of the method if local details are important. However, it might also be noted that 
the example cited contains a slope discontinuity. For cases with continuous slopes, this 
smoothing should be less severe. 

Low-Reduced- Frequency Aerodynamic Forces 

Comparisons with othe r analytical results and experiment.- The exact solutions for 
the low-frequency or stability-derivative -type aerodynamics of references 27 and 36 are 
available for a symmetrical wedge pitching about an axis located on the chord at any 
chordwise location. Direct algebraic comparisons with these results would be rather 
lengthy. The equations of reference 27 were programed and calculations made for an 
extensive range of M, y, and 0. The results obtained were identical to the results 
from using equations (14) and the transformation equations of appendix B. Furthermore, 
the results presented in the figures of reference 36 were also reproduced by using the 
one-strip equations. However, no results for Lj are available for comparison. 

Although not rigorously demonstrated, the results of the one-strip approximation (eqs. (14)) 
are thought to be exact for the low-frequency aerodynamics. Further comparisons with 
the linear theory of reference 14 and the second-order theory of reference 44 were also 
made by expanding the one-strip results in powers of © and obtaining identical results. 
Similar expansions for hypersonic small-disturbance conditions agreed with the results 
of reference 29. 

With the previously noted smoothing of the wave behavior by the method of integral 
relations, the degree of validity of the present results may seem surprising. It was found 
in reference 36, however, that linear functions in x and y satisfied the perturbation 
equations for the pitching wedge and the caret wing. The linear y -function is the form of 
solution assumed by the one-strip method of integral relations (eq. (6)) and the linear 
x-function is the form assumed for the low-frequency approximation (eq. (12)). 

Results calculated by using equations (14) are compared with the experimental data 
for a symmetrical pitching wedge (ref. 42) in figures 3 and 4. Generally good agreement 
is obtained, with the best agreement being for the thinner wedge and higher Mach number. 
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Selected analytical results .- The six normal-force coefficients (eqs. (13) and (14)) for 
a single inclined surface are presented in figures 5 and 6 for M = 2 and M = °°, respec- 
tively, and for several values of the ratio of specific heats y. All coefficients vary con- 
siderably with inclination angle, particularly for M = and approach infinity as detach- 
ment is approached. The coefficients and kL^ also change sign as detachment is 

approached. As the damping-in-pitch (about the leading edge) is related to kL^ (e.g., 
see eqs. (13)), the damping- in- pitch thus changes to an unstable value near detachment. 

This unstable pitch-damping has also been discussed in reference 36. 

The effect of y (figs. 5 and 6) is primarily to shift the detachment angle, and thus 
even small changes in the value of y have large effects near detachment. 

The fore and aft coefficients and kLg are, of course, zero for 0 = 0 and 
are somewhat small until the detachment angle is approached. 

It appears that the detachment condition leads to singular behavior in each stability 
derivative. That such is the case for the static forces is apparent from the charts of ref- 
erence 45 which show that dp g yci0 is infinite at detachment. For sonic local embedded 
flow the static derivatives are not infinite, but are extremely large. Furthermore, as dis- 
cussed in reference 46, the quantity D (see eq. (10)) in the denominator of the dynamic 
derivatives (see eqs. (14)) vanishes at the Crocco point, which occurs for a value of 0 
slightly larger than for sonic local flow and less than 0^. For the small range of 0 
for subsonic but attached flow, theories for an infinite wedge are no longer valid, as the 
subsonic flow field between the shock and the surface is not uniform for a finite wedge 
(ref. 47). Some implications of the approach to detachment on the flutter characteristics 
of diamond airfoils have been discussed in reference 25. 

Generally, in linear aerodynamics, the presence of a singularity such as that shown 
herein indicates a need for a more complete treatment of the nonlinear governing flow 
equations in order to describe a more realistic behavior. In the present paper, however, 
at least up to conditions of sonic local flow, accurate infinitesimal amplitude results are 
obtained. Furthermore, the usual transonic refinement to linear theory retains a local 
x-derivative of the steady flow field as a variable coefficient in the unsteady equations 
(ref. 48); herein, the term is zero. The present results suggest that infinitesimal motion 
has limited practical significance near detachment since even small finite amplitudes 
result in rapidly varying aerodynamic forces. Near detachment, the condition for infini- 
tesimal motion is that the amplitude of motion must be small compared with 0^ - 0, 
where 0^ is the detachment angle. 

Variation of Aerodynamic Forces With Reduced Frequency 

Comparisons with other results .- There are no results known to the author for the 
variation of the aerodynamic forces with reduced frequency as computed from the super- 
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sonic flow solution for the wedge of reference 26 or the extension of reference 27. Fur- 
thermore, such computations would be a rather lengthy task. Consequently, typical results 
of the integration of equations (7) are compared for 0 = 0 with supersonic linear theory 
for the flat plate (ref. 14), are compared for small 0 with second-order theory (ref. 44), 
and are compared with hypersonic small-disturbance theory (refs. 29 and 30) in figures 7 
to 9. In general, the one-strip integral method predicts the correct trends, but is accu- 
rate only when the frequency effects are small such as for thin surfaces at high Mach 
numbers (fig. 9) . Such a result might be anticipated in view of the previously discussed 
treatment of wave behavior by the integral method. 

Selected results .- The results of integrating equations (7) for M = y = 7/5, 
and 0 = 25° are given in figure 10. The frequency effects are small for k < 1; essen- 
tially no frequency effects are apparent for the usual range of interest for hypersonic 
speeds of k < 0.1. 

As previously discussed, the low-frequency aerodynamic coefficients vary rapidly 
with inclination angle as detachment is approached. Corresponding frequency effects are 
presented in figure 11. A very strong frequency dependence even at low reduced frequen- 
cies is evident as detachment is approached. In view of the previous comparisons with 
other theories (figs. 7 to 9), these results may be considered only qualitative in nature. 

The trends indicated are important ones, however. For example, the coefficient kL^ 
for 0 = 22° (fig. 11) varies so rapidly as to suggest that the linear results may be phys- 
ically unrealistic. If a nonlinear analysis were carried out for this region, then the phe- 
nomena of subharmonics or higher harmonics may occur which are rejected by the linear 
theory. In fact, the quasi-static results of references 34 and 39 indicate that for some 
conditions the unsteady pressure waveform is not a single harmonic of the frequency of 
oscillation and that shifts in the mean pressure level can occur even for moderate ampli- 
tudes of motion. Taking such nonlinear phenomena into account in a full dynamic analysis 
would require considerably more effort than would be required for a conventional linear 
analysis. 

Rapid variation of the coefficients with k may cause rapid variation of flutter 
characteristics with parameters that strongly influence the resulting k at flutter (such 
as the mass ratio; e.g., see ref. 49). 

CONCLUDING REMARKS 

A one-strip approximation of the method of integral relations has been applied to the 
perturbation equations describing small motions of an inclined flat surface. Accurate 
results were obtained with the method for the aerodynamic forces at low reduced fre- 
quency. Good agreement with available experimental data for a pitching wedge was also 
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obtained. The method, in the one-strip approximation at least, is of limited accuracy for 
cases in which the aerodynamic forces vary rapidly with reduced frequency. 

The results indicate that for the range of inclination angles near shock detachment, 
the aerodynamic forces vary rapidly both with inclination angle and with reduced frequency. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., January 26, 1972. 



APPENDIX A 


PERTURBATIONS OF FLOW VARIABLES AT A MOVING 
OBLIQUE SHOCK WAVE 


The incremental values Pp Pp Up and resulting from perturbations of 
slope and normal velocity of the moving oblique shock wave are required as boundary 
conditions for a perturbation analysis. The shock- wave perturbation is measured nor- 
mal to the steady shock wave, as shown in figure 1. The appropriate perturbations of 
the Rankine-Hugoniot conditions in a coordinate system alined with the steady shock wave 
have been developed in reference 26. They have also been given in the surface coordinate 
system used here (fig. 1) specialized to y = 7/5 (refs. 26 and 27), and for general val- 
ues of y but in somewhat different nondimensional units (refs. 26 and 36). Here, the 
shock derivatives are given for general y in the nondimensional units used for equa- 
tions (2) and are presented in much simpler algebraic form than those of references 26, 
27, and 36. 


With the complex amplitude of the normal velocity of the shock wave given by 


v n = iS^x) 


(Al) 


the complex amplitudes of the perturbations of the flow variables are written in the fol- 
lowing form: 

p 6 = V 1 + iP v 6 l 

p 6 = R^j + ilVj 

> (A2) 

u 6 = U ^1 + iU v 5 l 
v 6 = V^i + iV v 6j 

J 

where P^ = 8p g / 9 |3p P v = 9pg/8v n , and so forth, are evaluated at the steady shock wave 
position and subscript 6 denotes perturbation of conditions just aft of the shock wave. 

By letting 

M n = M oo sin % 1 


U 


n " 


2(M n 2 + 1) 
(y + l)M n 2 




(A3) 
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APPENDIX A — Concluded 


the shock derivatives can be written as 


4p 

p v = ; ^prr-sin/So 

^ + 1)p 0 u 0 


Rv = 


4p r 


(y + l)p 0O M n 2 sin /3 q 


U v = -U n sin X Q 


V v = U n cos X 0 


(A4) 


and 


p _ V °o P V c 

P /3 --=r- cos /3 0 


cos ^0 


Uo = n sin 0 - =P- P v cos X n 
P u n P v 0 

0 OO 


v TO u n 


PQ 


Vo = - °°" — cos © - =*- P v sin X 
Z 3 u A P v 




(A 5) 


These relations (eqs. (A4) and (A5)) reduce to those given in references 29 and 30 in the 
hypersonic small-disturbance limit. It can also be shown with a similar development 
that perturbations of the shock wave tangential to the steady shock position result in 
second- order perturbations in the flow variables. 

The steady-flow values of pressure, density, and velocities within the shock layer 
are given in reference 45 as functions of y, M^, and /3 q . The corresponding relation 
between j8q and 0 is also given as a cubic in sin2j8^ with coefficients as functions 
of y and ©. These relations are solved numerically for the steady-flow-field 
parameters. 
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APPENDIX B 

CALCULATION OF AERODYNAMIC COEFFICIENTS FOR A WEDGE 
FROM SURFACE COEFFICIENTS 

General Comments 

The unsteady aerodynamic coefficients for a single isolated surface have been 
derived in this report for rigid-body pitch and translation perpendicular and parallel to 
the surface. The coefficients and k were based on the length of the surface c, and the 
pitch axis was located at the leading edge. The coefficients for a symmetrical wedge at 
an angle of attack are given in this appendix in terms of the surface coefficients. The 
pitch axis is assumed to be on the wedge midplane or chord. The derivation of these 
relations is briefly outlined and the results are summarized. 

Force and Motion Transfer 

The forces and moments for a symmetrical wedge (subscript w) are related to the 
previously given forces and moments for a single surface by 
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APPENDIX B - Continued 


The perturbation displacements are related by 


hy - ±hy >w cos 6^ - h x>w sin 5 ^ x 0) w( cos 
h x = h x ,w cos 6w ± h yjW sin 6 ^ - x 0jW (sin 6 W )0 ) 

= = 


(B3) 


where again the upper sign pertains to the upper surface. 

The force and moment coefficients are defined by equations (11) for a single sur- 
face. By defining the coefficients for the wedge in a manner similar to equations (11), 
the coefficients for the surface and the wedge can be related through the use of equa- 
tions (Bl) for the dimensional forces. Then substituting equations (B2) and (B3) into the 
expression for the surface coefficients, the wedge coefficients can be equated to the cor- 
responding combination of surface coefficients. 


Coefficients for Wedge 

The results for the coefficients of lift and moment for the wedge are 









(B4) 
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APPENDIX B - Continued 


and 


M l,w = 


( M l,u + M l,l) + ( M 7,u + M 7,z) tan Sw]/ cos26 w 




M 2,w = [( M 2,u + m 2 ,l) + ( M 8,u + M 8,z) tan ^v]/ c os 2 5 w 
M 3 ; w=( M 3,u + M 3, Z )/ cos46 w 

m 4,w = ( M 4,u + M 4,z|/ cos ^ 6 w 

M 7,w = j( M 7,u " m 7 ,l) ~ ( M l,u “ M l,^) tan 6 w]/ cos2 5w 

M8,w = [(Mi, u - M' 8>1 ) - (M2, u - M 2,z)tan ^ 


(B5) 


'cos ^6, 


The relationship for reduced frequency for the wedge is 
k = cos 6 W 


k = cos 6 W 


(B6) 


It may also be noted that for zero angle of attack, the coefficients for the upper and lower 
surfaces are equal and thus the coefficients w , kLg w , , and kMg w are zero. 

Such would not be the case for nonzero angles of attack, however. 

For pitch-axis locations other than x Q w = 0, the transfer relations are found from 
the same procedure to be the conventional ones (e.g., see ref. 14) and are 

L l,w = L l,w 


l 2,w = L 2, 


w 


L 3,w = L 3, w " 2x o,w L l,w 
L 4,w = L 4,w " 2x o,w L 2,w 
L 7,w = L 7,w 
L 8,w = L 8,w 


(B7) 
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APPENDIX B - Concluded 


and 


M l,w M l,w " 2x o,w L l,w 


M 2,w M 2,w “ 2x o,w L 2,w 


^3,w ^3,w “ 2x o ,w( M l, w + ^3,w " 2x o,w^ J l,'w^ 


M 4,w M 4 ,w " 2x o,w( M 2, 


,w + L 4 ,w 2x o,w L 2,w) 


(B8) 


M 7,w " M 7,w " 2x o,w L 7,w 


M 8,w M 8,w ’ 2x o,w L 8,w 




27 



APPENDIX C 


DESCRIPTION OF COMPUTER PROGRAM 
General Comments 

A brief description of a FORTRAN computer program for numerically integrating 
the complex system of equations (7) to obtain the aerodynamic coefficients as a function 
of k is presented. A library routine INT2A for integrating a system of real, first- 
order, differential equations is used to integrate equations (7) after expanding into real 
form. The pitch and two translational motions are treated separately. 

The program was written for use on the Control Data series 6000 computer systems 
using the Langley Research Center versions of the RUN compiler and SCOPE 3.0 oper- 
ating system. Approximately 14 significant figures were used in the computations. The 
RUN compiler permits the use of multiple arithmetic statement on one card when sepa- 
rated by the character $. It might also be noted that the quantity 17770000000000000000g 
signifies an indefinite (undefined) quantity. 


Input 

Each case consists of a single card (80 characters) of identification for labeling the 
printout only and a list of variables in a NAMELIST called INPDATA. The FORTRAN 
variables and their definitions are as follows: 


FORTRAN variable 

Definition 

XM 

Moo 

THWD 

0, degrees 

G 

y 

XOK 

starting value of k for beginning 
numerical integration 

XSTOPK 

stopping value of k 

CIK 

k-increment for numerically integrating 
the differential equation 

SPECK 

k-increment for printing results 


The input cards for the results presented in figure 10 are listed below: 

******************************* COLUMN NUMBER ********************************** 

00000000011111111112222222222333333333344444444445555555555666666666677777777778 

12345678901234567890123456789012345678901234567890123456789012345678901234567890 

FREQUENCY EFFECTS AT M=1000.„ GAMMA=1.4„ SURFACE I NCL I NATI ON=25 DEG. 

$ INPDATA G=1.4, / XM=1000. # .THWD=25. / KOK=.0,01,XSTOPK=l. / CIK-.0Dl # SPECK=.01$ 
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APPENDIX C - Continued 
Output 

The program lists each aerodynamic coefficient for each value of k requested 
(by SPECK) and the real and imaginary parts of the perturbation variables. The coeffi- 
cients are also written on TAPE 7 in coded form for subsequent use by plotting programs. 


Listing of Program 


OVERLAY IPF1STRP, 0,0 I 

PROGRAM PF1STRP! INPUT=1, OUTPUT ,TAPE7, TAPE5= INPUT) 

********* * *************** ************************ ********** *********** *** 


* PROGRAM PF1STRP CALCULATES THE UNSTEADY FORCES ON AN OSCILLATING, * 

* INCLINED FLAT SURFACE IN A SUPERSONIC PERFECT GAS FLOW. * 

* THE LINEARIZED PERTURBATIONS ABOUT THE MEAN STEADY FLOW ARE * 

* CONSIDERED, WITH SOLUTION GIVEN BY A ONE-STRIP INTEGRAL METHOD. * 

* TWICE THE FORCES AND MOMENTS FOR A SINGLE SURFACE ARE PRINTED. * 


************************************************************************* 


C 


C 

c 

c 


COMMCN/BLOCK1 /CUVAR (9) , DER(9)*VAR(9) , PDB, RODB , UDB, VDB , XM2S , XBOS, 
* RB0S,CVSR,CPLGI,CPTHI,CTNL2,DVSDXR,DVSDXI, RVD8, PDV , UDV , VDV , RC NL 
+,DRO,DIO 

DIMENSION IDFNT18) ,DATE!2I , ELE1 { 8 ) , ELE2 ( 8 > , EPRVAL (9 ) 

EXTERNAL DERSUB.CHSUB 

NAMELlST/INPCATA/XM,THWD,G,XCK,XSTOPK ,C IK, SPECK 


110 F OR MAT ( 5EI6.81 

109 FORMAT ( IHI //* RIGID BODY FORE AND AFT TRANSLATION*) 

108 FORMAT! Z7X*K*AX*KSQL3P*6X*KLAP*4X*KSQM3P*6X*KMAP*9X*BR*8X*BI* 

* 8X*PR*3X*P I *8X*UR*8X*U I*8X*DR*9X*DI*/ ) 

107 FORMAT! 1HI//* RIGID BODY PITCH*) 

106 FORMAT! /7X*K*7X*L1*7X*KL2*7X*MIP*6X*KM2P*8X*BR*3X*BI* 

+ 8X*PR*3X*PI*8X*UR*8X*UI*8X*0R*8X*DI*/) 

108 FORMAT! 1H1//* RIGID BODY PLUNGE NORMAL TO SURFACE*) 

104 FORMAT! X F6 .4 , 4GI 1 .4 , 8G1 0 . 3 ) 

103 FORMAT! /7X*K*7X*L7*7X*KL 8*7X*M7P*6X*KM9P*8X*BR*8X*BI* 

+ 8X*PR*8X*P I*8X*UR*8X*UI *8 X*DR* 3X*D I */ ) 

102 FORMAT!//* PDB=*G16.8,* R0DB=*G16.8, * UDB=*G16.8,* VDB=*G16.8 
♦/* PDV=*G16 . 8 , * RODV=*G16.8,* UDV=*G16.8,* VDV=*G16.8/) 

101 FORMAT! /2X*M2=*G16. 8,* P21=*G18.8,* R21=*G16.8,* V21=*G16.8,* B0=* 
+G16.8) 

DATA FMT1,FMT2,FMT3/6H( 8A10) ,10H< 1H110A10) , 1CH( * I0B*I4)/ 

VAR! 1, ..,9)=X,BR,BI,PR,PI,UR,UI ,DR,DI. RF=.5*V21*X 


REWIND 7 

1 READ FMT1, IDENT $ I F ( EOF ,5 ) 999 , 2 

2 CALL DAYTIM!DATE) % PRINT FMT2 , I DENT , DATE $ READ INPDATA 
PRINT INPDATA $ THW=THWD/57. 2957795 130823 

CALL WEDGE! XM,Gt THW , XM2 , P21 , R21 ,V21 , BO, IEOBS) 

IF! IEOBS.EO.O)GO TO 3 S PRINT FMT3, IEOBS S GO TO 1 

3 PRINT 101»XM2,P21»R21»V21»B0 

CALL SHKDERV ! XM, G»THW , BO »R21 ,V2 1 , PDV, RODV, UDV , VDV , P DB , RODB , UDB , 
* VDB) $ XM2S=XM2**2 $ XBOS=XM2S-l. 

PRINT 102, PDB, RODB, UDB, VDB, PDV, RODV, UDV, VDV $ PRINT 105 
R BOS= 1 . /XBOS $ RVDB=1./VDB $ PC0N=-R21*V21 % C0NM=PC0N/3. 
CONTM=— 2.+R21/3. $ XO=2 . *XOK /V2 1 $ XSTOP=2 . *XS TOPK/ V2 1 
CI=2.*CIK/V2l $ SPEC=2.*SPECK/V21 
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APPENDIX C - Continued 


RIGID BODY PLUNGE 

CVSR=0. $ CPLGI=-1. % CPTH E=0. * 11=0 $ DRO=-COS ( BO-THW) $ D10=0. 
CALL INIC(DPDXR»DPDXI»XO,THW»BO» PROO, PIOO» TNL ) $ CTNL2=2./TNL 
CL 1Q=— R2 1*D PDXR $ CMlP0=4.*CL10/3. $ CKL20=PCON*PIOO $ XK=0. 
WRITE<7,li0)XK,CL10,CKL20,CMlP0,CKL20 
PRINT 106 SPRINT 104 *XK,CL10»CKL20»CM1PO»CKL20 

CALL INT2AI 1 1 , 8, DUM ,C I , SPEC , OUM ,DUM , VAR , CUV AR , DER , ELE l , ELE2 , DUM, 

+ ERRVAL, DERSUB, CHSUB, OUM) i XOO=VAR(1I 
C L 1=— R21* ( PR OO+VAR l 4 ) )/X00 $ CKL2= . 5*PC0N* ( PI OO+VAR ( 5 )) 

CM1 P=CONT M* ( PRGO+2 • *V AR ( 4 ) ) / XOO $ CKM2P=C0NM* ( PI 00+2 . *VAR ( 5 ) ) 

X SA V= VAR C 1 ) $ PR SAV=VAR ( 4 ) $ PISAV=VAR(5) $ XK= . 5 *V2 l*VAR ( l ) 

WRITE (7, 110) XK,CL1,CKL2,CM1P,CKM2P 

PRINT 104,XK,CL1,CKL2,CM1P,CKM2P, (VAR( I), 1=2, 9) $ 11=1 

7 CALL INT2AI 1 1 , 8 , DUM , C I , SPEC, DUM , DUM , VAR , CUV AR , DER ,ELE 1 , ELE2 , DUM, 

+ EPRVAL *DERSUB,CHSUR*DUM) 

X S= VAR ( 1 ) $ RXSQ=1. /XS**2 $ DX=XS-XSAV 
R X= XSAV/XS S RX2=RX*RX $ RX3=RX2*RX 
CL 1 =RX2*CL1-R21*DX*(PRSAV+VAR<4) )*RXSQ 
CKL2=RX*CKL2+.5*PC0N*0X*(PIS AV+VARI 5) >/XS 
TX1=VAR( l)+2.*XSAV $ TX2=2 .*VAR ( 1I+XSAV 
CKM2P=RX2*CKM2P+CCNM*RXSQ*DX*(TX1*P ISAV+TX2*VAR( 5 ) ) 
CM1P=RX3*CM1P+C0NTM*RXSQ*DX*(TX1*PRSAV+TX2*VAR(4) )/XS 
X SA V= VA R ( 1 ) $ PR SAV=VAR (4 ) $ PISAV=VAR(5) S XK=.5*V2l*VAR(l) 

WRITE (7, 110) XK,CL1,CKL2,CM1P,CKM2P 

PRINT 104,XK,CL1,CKL2,CM1P,CKM2P, <VAR( I), 1=2,9) 

I F( XS.LT .XSTOP)GO TO 7 $ ENDFILE 7 $ PRINT 1CP 

RIGID BODY FORE AND AFT TRANSLATION 

CVSP=CPLGI=CPTHI=DIO=0. $ 11=0 $ DRO= S I N ( BO-THW ) 

CALL INIC(DPDXR,DPDXI,XO,THW,BO,PROO,PIOO,TNL) $ PRINT 103 
CL70=-R21*0PDXR S C M7 P0=4 . +CL70/3 . $ CKL80= PC ON*P 100 $ XK=0. 

WRITE (7, 1 10 )XK,CL70,CKL80,CM7P0,CKL80 
PRINT 1P4,XK,CL70,CKL80,CM7PC,CKL80 

CALL INT 2 A ( I 1,8, DUM ,C I , SPEC , DUM , DUM , VAR , CUV AR , DER , ELE 1 , ELE 2 , DUM, 

+ EPRVAL, DER SUB, CHSUB, DUM) $ X00=VAR(1) 

CL7=-R21* ( PROO+VAR (4) )/X00 $ CKL8=. 5*PC0N*< PI 00+V AR ( 5 ) ) 
CM7P=C0NTM* ( PRCC+2.*VAR(4) )/X00$ CKM8 P=CONM*( P 1 00+2 . * VAR ( 5 i ) 

X SA V= VAR ( 1 ) $ PR SAV= VAR ( 4 ) $ PISAV=VAR(5) $ XK=. 5*V21*VAR ( 1 ) 

WRITE (7, 110) XK,CL7,CKL8,CM7P,CKM8P 

PRINT 104, XK , CL7 , CKL8 ,CM7P ,CKM3P, (VAR( I ( ,1=2,9) t 11=1 

8 CALL INT 2A ( 1 1 , 8 , DUM , C I , SPEC , DUM , DUM , V AR , CUV AR , DER , EL E 1 , ELE2 , DUM, 

+ ERRVAL,DERSUB, CHSUB, DUM) 

X S = VAR ( l ) S R XSQ= l • / XS**2 * OX = XS-XSAV 
R X= XSAV/XS t RX2=RX*R X $ RX3=RX2*RX 
CL7=RX2*CL7-R21*DX*(PRSAV+VAR(4) )*RXSQ 
CKL8 = RX*CKL8+.5*PC0N*DX*(PISAV+VAR( 5) )/XS 
TX1=VAR( 1)+2.*XSAV $ TX2=2.*VAR ( 1 )+XSAV 
CKM8P=PX2*CKM8P+CCNM*RXSQ*DX*(TX1*PISAV+TX2*VAR ( 5 ) ) 
CM7P=RX3*CM7P+C0NTM*RXSQ*DX*(TX1*PRSAV+TX2*VAR (4) ) /XS 
X SA V= VA R ( 1 ) $ PR SA V= VAR ( 4 ) $ PISAV=VAR(5) S XK=.5*V21*VAR(1) 

WRITE (7, LID )XK,CL7,CKL8,CM7P,CKM8P 

PRINT 104,XK,CL7,CKL8,CM7P,CKM8P,(VAR( I) ,1=2,9) 

I F ( XS.LT. XSTOP)GO TO 8 $ ENDFILE 7 $ PRINT 107 
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RIGID BODY PITCH 

C VS R=— 1 • $ CPLGI=0. $ CPTHI=-1. $ 11 = 0 $ DRO=DIO=C. 

CALL INIC(DPOXR,DPDXI f XO,THW,BO,PROO,PIOO,TNL)$ PRINT 108 
CKSL3PO=V21*PCON*PROO $CKL4P0=PCGN*DPDX I SC KM4P0=4 • *C KL4 PO/ 3 • 
XK= 0 • S WRITE(7,110)XK»CKSL3FQfCKL4P0»CKSL3P0,CKM4P0 
PRINT 104, XK,CKSL3P0,CKL4P0,CKSL3P0,CKM4P0 

CALL INT 2 A ( 11,8, DUM , C I , SPE C , DUM , DUM , V AR , CUV AR , DE R , E LE 1 , EL E2 , DUM 
+ ERRVAL,DERSUB»CHSUB,DUM) S X0C=VARC1) 

CKSL3P=* 5*PCON*V21* (PR0C+VARC4) I $ CKL4P=PC ON* ( P 1 00 + V AR ( 5 ) ) / XOO 
CKSM3 P=CQNM*V21* C PRO0+2 • *V AR ( 4 ) ) 

CKM4P=2**C0NM*( PI 00+2 • * V AR <5)1 /XOO 

X SA V= VA R C 1 I $ PR SAV = V AR ( 4 ) $ PISAV=VAR(5) S XK= . 5*V2 1 *VAR (1) 
WRITE(7 ,110) XK,CKSL3P,CKL4P,CKSM3P, CKM4P 

PRINT L04,XK,CKSL3P,CKL4P,CKSM3P,CKM4P, (VAR( I ) , 1=2,9) S 11=1 
9 CALL INT2A < I I , 8 , DUM , C I , SPEC , DUM , DUM , V AR , CUV AR , DE R , EL E 1 , EL E 2 , DUM 
+ ER RVAL , DER SUB ,CHSUB , DUM ) 

X S= VAR { 1 ) S R XSQ=1 • / X S**2 S CX=XS-XSAV 
RX=XSAV/XS S R X2 = RX*R X S RX3=RX2*RX 
CKSL3P=RX*CKSL3P+.5*DX*V2l*PCON*(PRSAV+VAR (4) ) /XS 
CKL4P=RX2*CKL4P +D X*R X SQ*PCO N* ( P I S AV + VAR ( 5 I ) 

TX1=VAR ( 1) +2.*XSAV S TX2 = 2 .*VAR (1) +XSAV 

CKSM3P=RX2*CKSM3P+C0NM*V21*RXSQ*DX*{TXl*PRSAV-»-TX2*VAR(M ) 
CKM4P=RX3*CKM4P +2.*C0N^*RXSQ*DX*(TX1*PISAV+TX2*VAR (5 ) )/XS 
X SA V= VA R ( 1 ) $ PRSAV = VAR(M $ PISAV=VAR<5> S XK= . 5*V2 1 *VAR ( 1 ) 

PRINT 104,XK,CKSL3P,CKL4P,CKSM3P,CKM4P, ( VAR ( I) , 1=2,9) 

WRITE(7 , 110 )XK,CKSL3P,CKL4P,CKSM3P, CKM4P 
IF( XS.LT.XSTOPIGO TO 9 S ENDFILE 7 $ GO TO 1 
999 REWIND 7 
END 


PROGRAM PF1STRP 


APPENDIX C - Continued 


SUBROUTINE INICf DPDXR , DPDX I , XO , THW, BO, PRO, PID, TNL ) 
************************************************************************* 

* SUBROUTINE INIC CALCULATES THE INITIAL VALUES ANO DERIVATIVES FOR * 

* BEGINNING THE NUMERICAL INTEGRATION AT XO. * 

************************************************************************* 

COMMON/BLOCK1/CUVARI9) ,DER!9) ,VAR!9) ,PDB,RODB,UDB,VDB,XM2S,XBOS, 

+ RBOS*CV$R,CPLGI»CPTHI»CTNL2»DVSDXR»DVSDXI*RVDB»PDV»UDV,VDV*RCNL 
+, DR 0,010 

VSR=CVSR $ VSI=CPLGI $ DVSDXI=CPTHI $ DVSDXR=0. 

TNL=TAN( BO-THWI * SNL=S IN ( BC-THW ) $ CNL=COS( BO-THW) $ RCNL=1./CNL 
BRO=(VSR+VOV*OIO)/VDB $ PRO=PDB*BRO-PDV*DIO 
BIO=( VSI-VDV*DRO)/VDB t PIO=PDB*BIO+PDV*DRO 
URO=UDB*BRO-UDV*DIO $ UIO=UDB*B IO+UDV*DRO 
PRINT 1 »VSRtVSI ,BRO,BIO,PRO, PI 0 ,URO ,UI 0, DRO , D I 0 

1 FORMAT!/* VSR=*G16.8,* VSI=*G16.8,* BR0=*G16.8,* BI0=*G16.8/ 

♦ * PR0=*GI6.8,* P 1 0=*G16 • 8 » * UR0=*G16.9,* UI0=*G16.8/, 

+ * DR0=*G16. 8 , * DI0=*G16.8/> 

BNUM1=1 «-XBOS*TKL*TNL % BNUM2=. 5*XB0S*TNL*TNL 
BNUM3=TNL*( XM2S*PDV-UDV+ . 5*V0V*XBOS*TNL ) 

BNUM4=TNL*( XM2S*PCB— UDB+VDV/ SNL+ ,5*XB0S*TNL*( VDB+2.*PDV/SNL) ) 
DEN0M=VDB*X80S*TNL*PDB 

DBDXR=( BNUM1*DVSDXR+BNUM2*VS I+BNUM3*DR0+BNUM4*BI0) /DENOM 

DBDX I=! BNUM1*0VSDX I— BNUM2*VS R+BNUM3*D 1 0— BNUM4*BR0 ) /DENOM 

DPDXR=PDB*DBDXR+TNL*!DVSDXR-.5*{ VS I +VDV*DRO+B I 0* I VDB+2 .*PDV/SNL > > I 

DPDXI=PDB*D8DXI+TNL*I DVSDXI+.5*! VSR-VDV*DIO+BRO*I VDB+2. *PDV/SNU ) I 

DUDXR=-OPDXR+UIO $ DUDX I -- DP CX I-URO 

PRINT 2, DBOXR,DBOXI ,DPDXR,DPDXI ,DUDXR,DUDXI 

2 FORMAT!* 08DXR=*G16 . 8 , * D BOX I=*G 16 . 8, * DPDXR= *Gl 6 . 8 , * DPDXI=*GI6.8 

+ ,/ * DUDXR=*G16 . 8 » * DUDX I=*G16. 9/ ) 

VAR 12 ) = 8R0-*-X0*DBDXR $ VAR ( 3)=BI0+X0*DBDXI i VAR(1)=X0 

VAR m=PR0«-X0*DPDXR $ VAR ( 5) -P I0+X0*DPDXI 

VAR !6 )=URO+XO*DUOXR $ VAR !7)=UI0+X0*DUDXI 

VAR 1 8 l = DRO*XG*BRO*RCNL $ VAR !9 I=DIO+XO*BIO*RCNL $ RETURN 

END SUBROUTINE INIC 

SUBROUTINE CHSUB $ RETURN 

******* ** ************* ************ *************************************** 

* DUMMY SUBROUTINE CALLED BY I NT2 A. * 

************************************************************************* 

END SUBROUTINE CHSUB 
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SUBROUTINE DERSUB 

***************************** ***** ********* *********** **************** *** 

* SUBROUTINE DERSUB CALCULATES THE FIRST ORDER DERIVATIVES FOR EACH 

* STEP OF THE NUMERICAL INTERGR AT I ON . DERSUB IS CALLED BY INT2A. 
************************************************************************* 

C0MM0N/BL0CK1/CUVAR(9),DER(9),VAR(9) , PDB, RODB , UDB , VDB , XM2S , XBOS, 

+ RBOS,CVSR,CPLGI,CPTHI,CTNL2,DVSDXR,DVSDXI , RVDB, PDV,UDV,VDV,RCNL 
+,DRO,DIO 

X= CUVAR (II $ BR=CUVAR<2) $ BI=CUVAR(3) $ PR=CUVAR<4) 

PI=CUVAR(5) $ UR=CUVAR ( 6 ) $ UI=CUVAR(7) $ DR=CUVAR(8> $DI=CUVAR<9) 

V SR=CVSR t VSI=(CPLGI+X*CPTHI ) t RX=I./X 
DER <21=1 (VDB*BR*CTNL2*(PR-PDB*BR)+(CTNL2*PDV-VDV)*DI-VSR)*RX 
l — D VSDXR+ ( ( VDB*VDV*RCNL)*BI+VSI+VDV*DR) )*PVDB 
DER <3 )= ( (VDB*BI+CTNL2*( PI-PDB*R I )- ( CTNL2*PD V- VDV ) *DR- VS I I *RX 
I -DVSDXI— ( ( VDB+VDV*RCNL )*BR+VSR— VDV*DI) )*RVDB 
DER(4)=(RX* ( CTNL2*VSR+ ( CTNL2*VDV— XBOS*PDV )*DI+(PDB*XBQS-VDB* 

1 CTNL2)*BR-XBOS*PR»-XBOS*PDB*DER< 2) +( ( PDV*XM2 S-UDV) *DR-UI 

2 +X M2S* PI + ( PDB*XM2S-UDB+XB0S*PDV*RCNL)*BI ) ) *RBOS 

DER (51= (RX*(CTNL2*VSI-(CTNL2*VDV-XB0S*PDV)*DR+<PDB*XB0S-V0B* 

1 CTNL2) *BI— XBOS*PI l-XBOS*PCB*DER ( 3) +( ( PDV+XM2 S-UDV) *DI+UR 

2 — X P2S*PR— ( PCR*XM2S-UDB*XBOS*POV*RCNL)*BR) )*RBOS 
DER (6) = RX*( C PDB+UDB ) * BR — PR— UR— (P0V+UDV)*DI )-DER(4) 

I -( PDB+UDB )*DER(2)+(UI+(UDB+(U0V+PDV)*RCNL)*BI+U0V*DR) 

DER (7 )=RX*( ( PCB+UDB ) * B I — P I — UI+( POV+UDV ) *DP ) -DE P ( 5 ) 

1 -(PDB+UDB)*DER(3)-<UR+(UDB-MU0V+PDV)*RCNL)*BR-UDV*DI ) 

DER (8 )=BR*RCNL t DE R ( 9 ) = B I *RCNL $ RETURN 
FND SUBROUTINE DERSUB 


SUBROUTINE SHKDF RV < XM ,G , THW, BO , ROO, UO , PV ,R OV, UV , VV, PB , ROB, UB , V B) 
************************************************************************* 

* CALCULATES SLOPE AND NORMAL VELOCITY DERIVATIVES OF P,RO,U, AND V FOR* 

* AN OBLIQUE SHOCK WAVE. DERIVATIVES ARE NORMALIZED BY FLOW VARIABLES* 

* BEHIND THE SHOCK WAVE. * 

********************************************************************** *** 


SLO=SIN( BO— THW) * C LO=C OS ( BO-THW ) 
SQMN= ( XM*SB ) **2 $ F Gl=4 . / ( G+l. ) 
PN= FGl*SB*RUO $ PV= PN/ROO 
ROV=FGl*ROO/( SQMN*SB) S 

UN=.5*FGI*( SQMN+1. ) /SQMN $ 

UV=— UN*SLO i VV=UN*CLO % 

RFTURN 

END SUBROUTINE SHKOERV 


t SB= SI N ( BO ) $ 

i RUO=l./UO 
$ PB=PV*CB*RUO 
ROB=ROV*CB 
UB=UN*RUO*SI N( THW ) 
VB=UN*RUO*CO S( THW) 


CB=COS( BQ) 


-PN*CLO 
-PN* SLO 
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SUBROUTINE WEDGE ( XM , G» THW, XM2, P20P1 , R2QR 1 ,V20Vl , BETAO, IER ) 
***************** ******************************************************** 

* CALCULATES SHOCK WAVE ANGLE BETAO (RADIANS) FOR A WEDGE OF ANGLE * 

* THETA (RADIANS) FOR INPUT VALUES OF MACH NUMBER XM AND GAMMA G - * 

* FLOW CONDITIONS BEHIND THE SHOCK ARE CALCULATED FROM BETAO, XM, AND G* 
************************************************************************* 

RINDF=OI777C00CO00OC0O00017 $ I F ( XM .GE . 1 . ) GO TO 1$ IER=1 $ GO TO 9 

1 IF(THW)8»7» 2 

2 STSQ=SIN(THW)**2 S SQM=XM*XM $ RSQM=1./SQM 

B=-RSQM*(SQM+2.)-G*STSQ $ D= ( STSQ-1 . ) *RSOM**2 

C=( 2.*SQM+1 . ) *PSCM**2+( .25*( G*( G+2. )+l . ) +RSQMM G- l. ) ) *STSQ 

CALL WDETACH(XM,G»TOET) $ IF ( THW.LE .TDET ) GO TO 3 t IER=3 $ GO TO 9 

3 X=-. 333333333333333*8 * TCORR=l. $ DO A 1=1,100 

C0RR=-(D+X*(C+X*(B+X)))/(C+X*(2.*B+3.*X)> S XT = X * X=X+CORR 
I F ( I .NE. I. AND. (TCORR*CORR.LE.O. .GR.X.EQ.XT) )G0 TO 5 

A TCORR=CORR $ IER=A t SMN=X*SQM $ GO TO 6 

5 IER=C $ SMN=X* SOM 

6 P20P1=(2.*G*SMN-G+1.)/(G+1.) $ R20R 1= ( G+ I . ) *SMN/ ( 2. + SMN* ( G-l . ) ) 

V20 V1=S0RT ( 1 A.*( SMN— 1 • ) * ( G+SMN+1 . ) / ( SMN*SQM* ( G*( G+2 . ) + 1 . > ) ) 
XM2=XM*V20V1*SCRT(R20R1/P20PL) $ BETAO= ASI N ( SORT ( X ) 1 S RETURN 

7 P20PI=R20R1=V20V1=1. $ XM2=XM $ 8ETA0=ASIN( I./XMI $ IER=0 $ RETURN 

8 I ER=2 

o P20P1=R20R1=V20V1=BETA0=XM2=RINDF * RETURN 
END SUBROUTINE WEDGE 

SUBROUTINE WDETACH( XM,G,DEL) $ XMS=XM*XM $ GP=G+1. 
************************************************************************* 

* CALCULATES WEDGE ANGLE FOR SHOCK DETACHMENT DEL (RADIANS) FOR INPUT * 

* VALUES OF MACH NUMBER XM AND GAMMA G * 

************************************************************************* 

XMNS=.25*(GP*XMS-A.+SQRT(GP*(XMS*(GP*XMS+8.*(G-1 . ) >+16.) ) >/G 
S INDEL=SQRT ( ( XMS-XMNS* ( 2 .*XMS+1 ,-XMNS* ( XMS + 2. -XMNS ) ) ) / 

+ ( l.+XMNS*(G-l. + .25*XMS*GP*GP-G*XMNS) ) )/XMi DEL = A SIN( S I NOEL ) $ RETU RN 
END SUBROUTINE WDETACH 
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SUBROUTINE INT2A( II , N , NT ,C I , SPEC , Cl MAX , IERR , VAR ,CUVAR, DER , ELE1 , 
IELE2 f ELT,EPRVAL,DERSUB,CHKSUR,ITEXT) 

A************************************************************************ 

* VERSION OF LRC CDC 6600 LIBRARY ROUTINE D2.4 * 

* USES FIXED INTERVAL, SINGLE PRECISION, 4TH-0RDER AOAMS-BA SHFORTH * 

* PREDICTOR, 4TH-0RDER ADAMS- MOLL TON CORRECTOR, 4TH-CR0ER RUNGE- 

* KUTTA STARTER. NT , C I M AX , IERR ,ELT , AND ITEXT ARE DUMMY (UNUSED). 

AA*A A* ************************************************* ********* ********* 

DIMENSION DER <21>,ELE1( 20) ,ELE2(20) ,ELT(20) ,ERRVAL(20) , TEMP (20) , 

1 DFRI(20),DER2(20),DER3(20),S1VAR(21),VAR(21) ,CUVAR(21) 

A=0.0 $ IF(II)1,I,2 

C INITIALIZATION SECTION 

1 IF(CI) 3, A, 3 

C SAVE Cl 

4 PRINT 1000 $ STOP 

1000 FORMAT! 11H0CI=0 STOP) 

3 H=C I 

IB I ERR=I S TO = VAR ( 1 ) t M0DE=1 $ 11= 1 $ N1=N+1 $ DO 5 J=l,NI 
CUV AR(J) = VAR(J) 

5 CONTINUE 

C EVALUATION SECTION HERE 

8 CALL DERSUB $ IF ( MODE .L E . 1 ) GO TO 7 S IF ( 1 1-3) 36,36,7 

36 CALL CHKSUB t IFUI.EQ.2) GO TO 1 

37 DO 38 J = I » N1 

38 VAP (J) = CUVAR( J) $ I F ( 1 1-3 ) 6, 7, 7 
7 RETURN 

6 I F ( SPEC ) 9,7,9 

9 DEL = VAR ( I ) -TO S DELP=DEL* ( 1 . + 1 . 0E- 1 1 ) 

I F ( ABS(DELP)— ABS(SPEC) ) 2,10,10 

10 TO = VAR ( 1 ) $ GO TO 7 

2 11=1 $ I F ( MCDE-4 ) 11,12,12 

C RIJNGE-KUTTA 

11 DO 20 J = 2 » N 1 $ DER3(J-1)=DER2( J-l) $ DER2 ( J-l ) =0ER1 ( J- 1 ) 

DER 1 ( J-l )=DER ( J ) t ELE1 ( J-l ) =DER ( J ) $ CUVAR(J)=A 

DELT=0.A*ELE1( J-l)*H $ SI VAR ( J ) =VAR ( J ) $ CUVAR ( J ) =S 1 VAR ( J ) +DELT 
20 CONTINUE $ SI VAR ( 1 ) =V AR ( 1 ) $ CUV AR ( 1 ) =S 1 VAR ( 1 ) +0 . 4*H % CALL DERSUB 
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I Ft 11-3)23,23,7 

23 CUVAP(1I=S1VAR(1) +0,4557 37 25421879*H $ DO 24 J=2,N1 
ELF 2 ( J- 1 ) = DFR ( J ) 

DEL T= ( 0 » 296 9776092477 5*ELE l(J-l) +0.15875 964 49 71 04*E LE 2 ( J- 1 ) ) *H 
CUVAR(J)=S1VAR(J)+DELT 

24 CONTINUE * CALL DERSUB * I F ( 1 1- 3 » 25,25,7 

2* CUVAR (I)=S1VAR(1)+H $ DO 26 J=2,N1 $ TEMPI J-l >=DER< J ) 

DELT=(0.21910C38822592*ELEi( J- U-3. 050965 14 86 9293*E LE 2 1 J- 1 ) 

1+3. 83286476 04 6701*TEMP( J-l) ) *H $ CUVAR ( J ) =S IV AR ( J ) +DE LT 

26 CONTINUE $ CALL DERSUB * I F< 1 1-3 ) 27, 27, 7 

27 DH=H $ CUVAR! 1>=VAR 1 1 )+DH $ DO 28 J=2,N1 

OOUB= I 0. 1747692 8226269*ELE II J-l >-0 . 551 48066287873*EL E 2 1 J-L ) 
1+1. 2055355°9396?2*TEMP( J-l )+C. 171184781 2 1952*DER I J) ) 

CUVAR (J ) =VAR(J)+DH* DOUB 

28 CONTINUE $ MODE=MODE+ 1 $ GO TO 8 

C ADAMS-MOULTON, ADAMS-BASHFORTH PREDICTOR 

12 CUVAR (1 >=VAR (l)+H $ DH=H/24.0 * DO 13 J=2,N1 

DOUB= (55.0*DER(J)-59.0*0ERl(J-l)+37.0*DER2( J- l I -9 . 0*DER3 I J-l ) ) 
CUVAR ( J)=VAR( J)+DH*DOUB 

13 CONTINUE $ DO 14 J=l,N t DER 3( J )=DER2(JI t DER2(J )=DER1(J) 

1* DER 1 ( J ) = DER ( J + l ) * CALL DERSUB * IF ( 1 1-3 ) IE , 1 5 , 7 

C ADA MS-MOULTCN CORRECTOR 

15 DO 16 J = 2 » N 1 $ TEMP=CUVAR ( J J 

DOU B= (9.0*DER( J >+l°.*DERK J-1)-5.0*DER2( J-1)+DER3( J-l) ) 

CUVAR (J )=VAR I J)+DH*DOUR 

16 FRRVAH J-l ) = I TEMP- CUVAR(J)) /14. 21052631579847 

19 GO TO 8 

FND SUBROUTINE INT2A 
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(a) t = 0.16; Mach 1.75. 


(b) r = 0.16; Mach 2.47. 


1- strip integral method 

O Experiment (ref. 42) 



(c) r = 0.24; Mach 1.75. (d) r = 0.24; Mach 2.47. 

Figure 3.- Variation of static moment derivative with pitch-axis location. 
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(a) t = 0.16; Mach 1.75. (b) r « 0.16; Mach 2.47. 



(c) r = 0.24; Mach 1.75. (d) r = 0.24; Mach 2.47. 

Figure 4.- Variation of pitch-rate-damping moment derivative with pitch-axis location 
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Figure 5.- Coefficients for a single surface at a Mach number of 2. k = 0. 
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1- strip integral method 

Second- order theory (ref. 44) 
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Figure 8.- Comparison of frequency effects for a wedge at 
a Mach number of 10/7. r = 0.10. 
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